#ifdef __MAKECINT__
#pragma link C++ class std::vector<int>+;
#pragma link C++ class std::vector<float>+;
#pragma link C++ class std::vector < std::vector<float> >+;
#pragma link C++ class std::vector<int>+;
#pragma link C++ class std::vector < std::vector<int> >+;
#pragma link C++ class std::string+;
#pragma link C++ class std::vector<string>+;
#endif

#ifndef FILTERS_H
#define FILTERS_H

#include <vector>
#include "utility.h"
#include "modifiers.h"
#include "jetclustering.h"
#include "observables.h"
#include <stdio.h>
#include <iostream>
#include "mt2_bisect.h"
#include "TRandom1.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TStopwatch.h"
#include "TF1.h"
#include "TUnuran.h"
#include "TUnuranContDist.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TFrame.h"
#include "TH1F.h"
#include "TBenchmark.h"
#include "TRandom.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TObject.h"
#include "services.h"

using namespace std;
using namespace mt2_bisect;

char blue[] = { 0x1b, '[', '1', ';', '3', '4', 'm', 0 };
char ublue[] = { 0x1b, '[', '4', ';', '3', '4', 'm', 0 };
char magenta[] = { 0x1b, '[', '1', ';', '3', '5', 'm', 0 };
char red[] = { 0x1b, '[', '1', ';', '3', '1', 'm', 0 };
char normal[] = { 0x1b, '[', '0', ';', '3', '9', 'm', 0 };

/**********************************************************/
// PURPOSE: A filter has as output a vector<vector<int> >
// which is the list of the collection of particles
// to be used to compute the observables
// each list might be a jet according to some jet algorithm
// or any other kind of cluster, as for instance 4 particles might
// be divided into two clusters with nearly equal mass or with similar deltaR
/**********************************************************/

typedef long double (*FUN)(vector<vector<int> >, vector<TLorentzVector>,
		vector<TRootLHEFParticle*>);

//**************************************************
extern "C" vector<vector<int> > filter_keeppartons(Int_t nmax,
		vector<TRootLHEFParticle*> partvector, int minstatus, int maxstatus) {
	//**************************************************
	vector < vector<int> > output;
	output.resize(0);
	vector<int> partons;
	partons.resize(0);

	for (int i = 1; i < 7; i++) {
		partons.push_back(i);
		partons.push_back(-i);
	}

	partons.push_back(21);

	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {
		int PID;
		PID = partvector[jcount]->PID; //each time changes
		if (isinclass(partons, PID)) {
			int status;
			status = partvector[jcount]->Status;
			if (isinrange(status, minstatus, maxstatus) == 1) {
				vector<int> temp;
				temp.resize(0);
				temp.push_back(jcount);

				output.push_back(temp);
			}
		}
		/*if ( !isinclass(classpid,PID)) {
		 cout<<"no interesting particle found"<<endl;;
		 }*/

	}

	return output;
}

//*************************************************
extern "C" vector<vector<int> > filter_trivial(vector<int> classpid,
		vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, int minstatus, int maxstatus,
		float mcut, TFile *f, TTree * t, int eventnum,
		vector<int> & reusedclustering, float alg) {
	//*************************************************
	Int_t nmax;

	//	TFile *f;
	//	f = filein;
	nmax = vec_partvector.size();

	vector < vector<int> > classmembers;
	classmembers.resize(0);

	int jetpid;
	jetpid = 93;

	// 94 95 96 97 are handled separately as they use the MCtruth

	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {
		int PID;
		PID = partvector[jcount]->PID; //each time changes
		if (isinclass(classpid, PID)) {

			//			if (isinclass(classpid, 0)) {
			//				cout << "il PID di " << jcount << " era " << PID << "ma l'ho preso per via dello zero" << endl;
			//			}

			int status;
			status = partvector[jcount]->Status;
			if (isinrange(status, minstatus, maxstatus) == 1) {
				vector<int> temp;
				temp.resize(0);
				temp.push_back(jcount);

				classmembers.push_back(temp);
			}
		}
		/*if ( !isinclass(classpid,PID)) {
		 cout<<"no interesting particle found"<<endl;;
		 }*/

	}

	vector < vector<int> > jadejets;

	/* ***************************** */
	/* CHECK IF IS NEEDED TO CLUSTER */
	/* ***************************** */
	if (isinclass(classpid, jetpid)) {

		/* faccio i jet */

		//		double mcut;
		//		mcut = 80 / 2.0; /* best ratio 2jet/3jet */
		//		mcut = 100.0;
		bool event_exists;
		event_exists = false;

		Int_t rootindex;

		//		if (file_exists == true) {

		//		cout << "/* check if the event has been already clustered */" << endl;

		event_exists = findeventnum(eventnum, rootindex, f, "tvec",
				"eventnumber");

		//		hello();

		//		cout << event_exists << " found event " << endl;
		//		}

		if (event_exists == true) {
			/*	*********************************** */
			/* READ FROM THE FILE */
			/*	*********************************** */
			// cout << "found event " << rootindex << " already clustered, gnam gnam" << endl;
			// puts("leggo i cluster dal rootfile");
			/*	*********************************** */
			/* THE TTree HAS BEEN PASSED AS POINTER */
			/*	*********************************** */
			//			 TTree *t;
			//			f->GetObject("tvec", t);
			/*	*********************************** */
			/* THE TTree HAS BEEN PASSED AS POINTER */
			/*	*********************************** */
			std::vector < std::vector<int> > *vpx = 0;
			TBranch *bvpx = 0;
			t->SetBranchAddress("clusters", &vpx, &bvpx);

			Long64_t tentry = t->LoadTree(rootindex);
			bvpx->GetEntry(tentry);

			jadejets = *vpx;
			reusedclustering.push_back(tentry);

			//printmatrix(jadejets);

			// Since we passed the address of a local variable we need
			// to remove it.
			t->ResetBranchAddresses();

			//			f->Close("R");

			//			delete f;

		} else {
			/* ***************************** */
			/* DO THE CLUSTERING             */
			/* ***************************** */
			vector<int> tobeclustered;
			/* ***************************** */
			/* ***************************** */
			/* TAKE CARE OF THE FIXED STATUS */
			tobeclustered = matrix2vector(
					filter_keeppartons(nmax, partvector, 1, 100));
			/* ***************************** */
			/* ***************************** */

			if (alg == 0) {

				jadejets = jade(vec_partvector, partvector, tobeclustered,
						mcut);
			} else if (alg == 1) {

				vector < TLorentzVector > modified_vectors;

				/* *************************************** */
				/*   PARAMETERS OF A PARTICULAR MODIFIER   */
				/* *************************************** */
				//				long double theta;
				//				theta = acos(-1.0);
				//				theta = acos(cosijMotherChildWPlus(Wclassmembers, vec_partvector, partvector));
				//				long double gamma;
				//				gamma = 1000.0 / 80.0;
				//				gamma = 2.5;
				/* *************************************** */
				/*   MODIFIERS OF THE PARTICLES 4VECTORS   */
				/* *************************************** */

				//				modified_vectors = modifier_boost(gamma, theta, partvector, vec_partvector);
				//				modified_vectors = modifier_boost_relative(gamma, theta, partvector, vec_partvector);
				//				modified_vectors = modifier_boost_identity(gamma, theta, partvector, vec_partvector);
				/* *************************************** */
				/*         END OF MODIFIERS                */
				/* *************************************** */

				modified_vectors = vec_partvector;

				jadejets = kT(modified_vectors, partvector, tobeclustered,
						mcut);

			}

			/* append this clustering to root file */
			/*	*********************************** */
			/* THE TTree HAS BEEN PASSED AS POINTER */
			/*	*********************************** */
			//			TTree *t = (TTree*) f->Get("tvec");
			/*	*********************************** */
			/* THE TTree HAS BEEN PASSED AS POINTER */
			/*	*********************************** */

			int *evn = 0;
			evn = &eventnum;

			cout << "event #" << eventnum << endl;

			t->SetBranchAddress("eventnumber", &eventnum); //, &bvpx);

			//TBranch *bvpx = 0;
			std::vector < std::vector<int> > *vpx = 0;
			vpx = &jadejets;

			t->SetBranchAddress("clusters", &vpx); //, &bvpx);

			t->Fill();

			//			t->Write(0, TObject::kWriteDelete, 0);
			// ipse dixit
			/* qui si lamenta */
			// 	Error in <TFile::WriteTObject>: Directory pairs_in_N.root is not writable
			//			t->Write(0, TObject::kWriteDelete, 0);
			//			t->Write();
			//			hello();
			//			cout<< "sto scrivendo il file"<<endl;
			//			f->Write();
			//			t->Write();
			//			t->ResetBranchAddresses();
			/* THIS WRITES TO THE FILE AFTER EACH CLUSTERING IS FINISHED */
			//			f->Write(0, TObject::kWriteDelete, 0);
			/* THIS WRITES TO THE FILE AFTER EACH CLUSTERING IS FINISHED */

			//			f->Close();
			//			delete f;
		}

	}

	/* ***************************************** */
	/* CHECK IF CLUSTERING ACCORDING TO MC TRUTH */
	/* ***************************************** */
	int mcjetpid;
	mcjetpid = 94;
	int mcpositivejetpid;
	mcpositivejetpid = 95;
	int mcnegativejetpid;
	mcnegativejetpid = 96;
	int mcbothjetpid;
	mcbothjetpid = 97;
	int mcpositivejetpidpy6;
	mcpositivejetpidpy6 = 99;
	int mcnegativejetpidpy6;
	mcnegativejetpidpy6 = 98;
	int mcbothjetpidpy6;
	mcbothjetpidpy6 = 100;

	if (vectorcontains(classpid, mcjetpid)
			|| vectorcontains(classpid, mcbothjetpid)
			|| vectorcontains(classpid, mcnegativejetpid)
			|| vectorcontains(classpid, mcpositivejetpid)
			|| vectorcontains(classpid, mcbothjetpidpy6)
			|| vectorcontains(classpid, mcnegativejetpidpy6)
			|| vectorcontains(classpid, mcpositivejetpidpy6)) {

		/* ******************************** */
		/* DO THE CLUSTERING USING MC TRUTH */
		/* ******************************** */
		vector<int> tobeclustered;
		/* ****************************************** */
		/* ****************************************** */
		/* TAKE CARE OF THE FIXED STATUS              */
		/* NOW IT RETURNS ALL THE FINAL STATE PARTONS */
		tobeclustered = matrix2vector(
				filter_keeppartons(nmax, partvector, 1, 100));
		/* ****************************************** */
		/* ****************************************** */

		int type;
		if (vectorcontains(classpid, mcjetpid)) {
			type = 94;
		}

		if (vectorcontains(classpid, mcpositivejetpid)) {
			type = 95;
		}
		if (vectorcontains(classpid, mcnegativejetpid)) {
			type = 96;
		}

		if (vectorcontains(classpid, mcbothjetpid)) {
			type = 97;
		}
		if (vectorcontains(classpid, mcpositivejetpidpy6)) {
			type = 99;
		}
		if (vectorcontains(classpid, mcnegativejetpidpy6)) {
			type = 98;
		}
		if (vectorcontains(classpid, mcbothjetpidpy6)) {
			type = 100;
		}
		//		jadejets = mcclusters(vec_partvector, partvector, tobeclustered, mcut);
		jadejets = mcclusters(vec_partvector, partvector, tobeclustered, mcut,
				type);

	}

	//printmatrix(jadejets);

	return mergevector(classmembers, jadejets);

}

//*************************************************
extern "C" vector<vector<int> > pick_list_in_OBS(vector<vector<int> > kplets,
		vector<vector<int> > classmemebers,
		vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, string parameters) {

	//*************************************************
	// the input string "parameters" must be a comma separated list without spaces
	// of the form variable1=val1,variabe2=val2
	// at run-time the use can specify an arbtrary number of variables
	// but each of them must be explictly declared in the code of the filter he/she is using
	// calling filter_pTmin(dummy=value)
	// assumes that the declaration of a variable "float dummy;" is present
	// at an appropriate point in the code of the filter.
	//*************************************************

	//*************************************************
	// this filter eliminate all the kplets but the one
	// that is the n-th in the sorting according to pT
	// pT can actually be any observable
	//*************************************************

	vector < vector<int> > outputkplets;
	outputkplets.resize(0);
	vector<float> object_to_sort;
	int maxrequested;

	//	float nth = get_float_parameter(parameters, "nth");
	string selection = get_string_parameter(parameters, "list");
	vector<int> selected_ns = get_separated_values(selection, ';');

	maxrequested = maxvector(selected_ns);

	//cout << "massimo jet richiesto: " << maxrequested << endl;
	//cout << "lista dei cosi richiesti" << endl;
	//printvector(selected_ns);

	FUN thisobs_func = get_the_obs_function(parameters);
	string OBS = get_string_parameter(parameters, "OBS");

//	string arranger = get_string_parameter(parameters, "arranger");
//	vector < string > arrangers;
//	arrangers.resize(0);
//	arrangers.push_back(arranger);

	// build the necessary pairs of pair<pT,EID> and sort according to the long double pT

	//cout << endl << "prima del sort" << endl;	printmatrix(previousclassmemeber);	cout << endl;

	object_to_sort = computeOBSonthekpletsimple(partvector, vec_partvector,
			kplets, classmemebers, thisobs_func);

//	printmatrix(kplets);

//	printvector(object_to_sort);	cout << endl;

	if (maxrequested <= object_to_sort.size()) {

		vector < vector<int> > sorted_kplets;
		sorted_kplets = sort2ndaccordingto1st(object_to_sort, kplets);

//		cout << "dopo il sort" << endl; printmatrix(sorted_kplets); cout << endl;

//		cout << "dopo il sort" << endl;	printmatrix(sorted_classmembers);	cout << endl;

		//	// print the sorted observables for a check
		//	for (Int_t jcount = 0; jcount < sorted_classmembers.size(); jcount++) {
		//
		//		vector<vector<int> > matrixtoadd;
		//		matrixtoadd.resize(0);
		//		matrixtoadd.push_back(sorted_classmembers[jcount]);
		//
		//		long double obs;
		//		// call the observable that is needed
		//		cout << pTofn(matrixtoadd, vec_partvector, partvector) << " ";
		//	}

		// strip all the particle except the one that satisfies the ordering requirement
		for (Int_t jcount = 0; jcount < sorted_kplets.size(); jcount++) {
			//		if (jcount + 1 == nth) {
			//			outputkplets.push_back(sorted_kplets[jcount]);
			//		}
			if (vectorcontains(selected_ns, jcount + 1)) {
				outputkplets.push_back(sorted_kplets[jcount]);
			}

		}

		return outputkplets;
	} else {
		//	cout << " too few objects to sort" << endl;
		outputkplets.resize(0);
		return outputkplets;

	}

}

//*************************************************
extern "C" vector<vector<int> > usePIDs(vector<vector<int> > kplets,
		vector<vector<int> > previousclassmemeber,
		vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, string parameters) {
	//*************************************************

	vector < vector<int> > outputclassmembers;
	outputclassmembers.resize(0);
	outputclassmembers = previousclassmemeber;

	vector<int> allowed_PIDs = get_CSV(parameters);

	//---->loop over the (clusters of) particles of the input list !!!from the end!!!<------
	for (Int_t jcount = previousclassmemeber.size() - 1; jcount > -1;
			--jcount) {
		for (Int_t kcount = previousclassmemeber[jcount].size() - 1;
				kcount > -1; --kcount) {

			//do the checks that the user wants to do to select particles

			int thisPID;

			thisPID = partvector[previousclassmemeber[jcount][kcount]]->PID; //each time changes
			// do the check on the computed value
			if (isinclass(thisPID, allowed_PIDs) == false) {
				outputclassmembers[jcount] = stripelement(kcount,
						outputclassmembers[jcount]);
			}
		}
	}

	return outputclassmembers;

}

//*************************************************
extern "C" vector<vector<int> > useTwoClassesOfPIDs(vector<vector<int> > kplets,
		vector<vector<int> > previousclassmemeber,
		vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, string parameters) {
	//*************************************************
// example:	useTwoClassesOfPIDs(1,2,3,4,-1,-2,-3,-4,21;5,-5)

	vector < string > PIDsClasses = get_separated_strings(parameters, ';');
	vector<int> allowed_PIDsOne = get_CSV(PIDsClasses[0]);
	vector<int> allowed_PIDsTwo = get_CSV(PIDsClasses[1]);

//	cout << "first sub-class:";
//	printvector(allowed_PIDsOne);
//	cout << "second sub-class:";
//	printvector(allowed_PIDsTwo);

	//cout << " filter by 2xPID " << endl;
	//printmatrix(previousclassmemeber);
	//printmatrix(kplets);

	vector < vector<int> > outputkplets;
	outputkplets.resize(0);
	outputkplets = kplets;

	for (Int_t kpletcount = kplets.size() - 1; kpletcount > -1; --kpletcount) {
//		scan from the end as not to modify the positions while removing stuff

		bool keepthiskplet = false;

		if (kplets[kpletcount].size() == 2) {

			int firstPID = partvector[previousclassmemeber[kplets[kpletcount][0]
					- 1][0]]->PID;

			int secondPID =
					partvector[previousclassmemeber[kplets[kpletcount][1] - 1][0]]->PID;

	//		cout << kplets[kpletcount][0] << " => "
		//			<< previousclassmemeber[kplets[kpletcount][0] - 1][0]
			//		<< " => " << firstPID << endl;

			//cout << kplets[kpletcount][1] << " => "
				//	<< previousclassmemeber[kplets[kpletcount][1] - 1][0]
					//<< " => " << secondPID << endl;

			if ((isinclass(firstPID, allowed_PIDsOne) == true)
					&& (isinclass(secondPID, allowed_PIDsTwo) == true)) {
				// the two are in differnt sub-classes

				//cout << "11" << endl;

				keepthiskplet = true;
			}

			if ((isinclass(firstPID, allowed_PIDsTwo) == true)
					&& (isinclass(secondPID, allowed_PIDsOne) == true)) {
				// the two are in differnt sub-classes
				keepthiskplet = true;

			//	cout << "22" << endl;

			}

			if (keepthiskplet == false) {
				outputkplets = stripelement(kpletcount, outputkplets);
			}

		} else {
			cout << " useTwoClassesOfPIDs works only for pairs " << endl;
		}

	}
	//---->loop over the (clusters of) particles of the input list !!!from the end!!!<------
//	for (Int_t jcount = previousclassmemeber.size() - 1; jcount > -1;
	//		--jcount) {
	//	for (Int_t kcount = previousclassmemeber[jcount].size() - 1;
	//		kcount > -1; --kcount) {

	//do the checks that the user wants to do to select particles

	//int thisPID;

	//thisPID = partvector[previousclassmemeber[jcount][kcount]]->PID; //each time changes
	// do the check on the computed value
	//if (isinclass(thisPID, allowed_PIDsOne) == false) {
	//		outputclassmembers[jcount] = stripelement(kcount,
	//			outputclassmembers[jcount]);
	//}
	//}
	//}

	return outputkplets;

}

//*************************************************
extern "C" vector<vector<int> > with_OBSminmaxrange(
		vector<vector<int> > inkplets, vector<vector<int> > inclassmemeber,
		vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, string parameters) {

	//*************************************************
	// the input string "parameters" must be a comma separated list without spaces
	// of the form variable1=val1,variabe2=val2
	// at run-time the use can specify an arbtrary number of variables
	// but each of them must be explictly declared in the code of the filter he/she is using
	// calling filter_pTmin(dummy=value)
	// assumes that the declaration of a variable "float dummy;" is present
	// at an appropriate point in the code of the filter.
	//*************************************************
	FUN thisobs_func = get_the_obs_function(parameters);
	string OBS = get_string_parameter(parameters, "OBS");
	string minmaxrange = get_string_parameter(parameters, "minmaxrange");

	vector<float> vresults;

	//	vector<string> assignments = getCSstrings(parameters);
	//	map<string, float> mymap;
	//	//	printvector(assignments);
	//
	//	mymap = assignment2map(assignments);
	//	// from here on the variables are initialized and valorized from the string provided as input

	float OBSmin, OBSmax;

	if (minmaxrange == "min") {
		OBSmin = get_float_parameter(parameters, "OBSmin");
		//				valorize(OBSmin, "OBSmin", mymap);
//		cout << OBS << " will be enforced to be larger than " << OBSmin << endl;

	} else if (minmaxrange == "max") {
		//		valorize(OBSmax, "OBSmax", mymap);
		OBSmax = get_float_parameter(parameters, "OBSmax");
//		cout << OBS << " will be enforced to be smaller than " << OBSmax
		//	<< endl;

	} else if (minmaxrange == "range") {
		OBSmin = get_float_parameter(parameters, "OBSmin");
		OBSmax = get_float_parameter(parameters, "OBSmax");

		//		valorize(OBSmin, "OBSmin", mymap);
		//		valorize(OBSmax, "OBSmax", mymap);
		//	cout << OBS << " will be enforced to be larger than " << OBSmin << endl;
		//cout << OBS << " will be enforced to be smaller than " << OBSmax
		//	<< endl;

	}

	//cout << "dim lista kplets " << inkplets.size() << endl;

	if (!inkplets.empty()) {

		vector < vector<int> > outputlist;
		outputlist.resize(0);

		outputlist = inkplets;
		//cout << "dim lista kplets " << inkplets.size() << endl;
		// if this is run as a filter for kplets it is exactly what you want
		// when run as filter on single particles it make sense anyways because in that case
		// the two input lists are equal

		// kplet mode
		//	vresults = computeOBSonthekpletsimple(partvector, vec_partvector, inkplets,
		//		inclassmemeber, OBS);

//	string arranger = get_string_parameter(parameters, "arranger");
//	vector < string > arrangers;
//	arrangers.resize(0);
//	arrangers.push_back(arranger);

		//	cout << "about to compute on "<< inkplets.size() << endl;

		vresults = computeOBSonthekpletsimple(partvector, vec_partvector,
				inkplets, inclassmemeber, thisobs_func);
//		cout << "done on "<< inkplets.size() << endl;

		for (Int_t jcount = inkplets.size() - 1; jcount > -1; jcount--) {
			//	cout << " checking the kpletlist's element " << jcount << endl;
			// do the check on the computed value
			//		if ( vresults[jcount] < OBSmin) {
			if (!meets_the_requirement(vresults[jcount], minmaxrange, OBSmin,
					OBSmax)) {
				//	cout << jcount << " must be removed because " << OBS << " is "
				//				<< vresults[jcount] << endl;
				//	cout << " the list of kplets before the removal is" << endl;
				//	printmatrix(outputlist);
				outputlist = stripelement(jcount, outputlist);
				//	cout << jcount << " has been removed because " << OBS << " was "
				//		<< vresults[jcount] << endl;
				//	cout << "now the list of kplets is" << endl;
				//	printmatrix(outputlist);
			}
		}
//		cout<<" done non empty"<<endl;
		return outputlist;
	} else {
		//	cout<<" done with empty"<<endl;
		return inkplets;

	}
}

//*************************************************
extern "C" vector<vector<int> > take_all(vector<vector<int> > inkplets,
		vector<vector<int> > previousclassmemeber,
		vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, string parameters) {
	//*************************************************

	if (inkplets == previousclassmemeber) {
		return previousclassmemeber;
	} else {
		return inkplets;
	}

}

////*************************************************
//extern "C" vector<vector<int> > with_minv_in_range(vector<vector<int> > kplets, vector<vector<int> > previousclassmemeber, vector<TLorentzVector> vec_partvector, vector<TRootLHEFParticle*> partvector, string parameters, FUN OBS) {
//
//	//*************************************************
//	// the input string "parameters" must be a comma separated list without spaces
//	// of the form variable1=val1,variabe2=val2
//	// at run-time the use can specify an arbtrary number of variables
//	// but each of them must be explictly declared in the code of the filter he/she is using
//	// calling filter_pTmin(dummy=value)
//	// assumes that the declaration of a variable "float dummy;" is present
//	// at an appropriate point in the code of the filter.
//	//*************************************************
//
//	vector<vector<int> > outputclassmembers;
//	outputclassmembers.resize(0);
//	outputclassmembers = previousclassmemeber;
//
//	vector<string> assignments = getCSstrings(parameters);
//	map<string, float> mymap;
//	mymap = assignment2map(assignments);
//	// from here on the variables are initialized and valorized from the string provided as input
//	float Mmin, Mmax;
//	valorize(Mmin, "Mmin", mymap);
//	valorize(Mmax, "Mmax", mymap);
//
//	//---->loop over the particles of the input list !!!from the end!!!<------
//	for (Int_t jcount = previousclassmemeber.size() - 1; jcount > -1; --jcount) {
//		cout << "******************* " << previousclassmemeber[jcount].size() << endl;
//		//do the checks that the user wants to do to select particles
//
//		vector<vector<int> > matrixtoadd;
//		matrixtoadd.resize(0);
//		matrixtoadd.push_back(previousclassmemeber[jcount]);
//
//		long double minv;
//
//		// call the observable that is needed
//		minv = minvofn(matrixtoadd, vec_partvector, partvector);
//		cout << "the mass of the pair ";
//		printmatrix(matrixtoadd);
//		cout << " was " << minv << endl;
//		// do the check on the computed value
//		if (minv < Mmin || minv > Mmax) {
//			outputclassmembers = stripelement(jcount, outputclassmembers);
//			cout << " this was *not kept* for the plot " << endl;
//		} else {
//			cout << " this was *kept* for the plot " << endl;
//		}
//
//	}
//
//	return outputclassmembers;
//
//}

////*************************************************
//extern "C" vector<vector<int> > with_etamin(vector<vector<int> > kplets, vector<vector<int> > previousclassmemeber, vector<TLorentzVector> vec_partvector, vector<TRootLHEFParticle*> partvector, string parameters, FUN OBS) {
//
//	//*************************************************
//	// the input string "parameters" must be a comma separated list without spaces
//	// of the form variable1=val1,variabe2=val2
//	// at run-time the use can specify an arbtrary number of variables
//	// but each of them must be explictly declared in the code of the filter he/she is using
//	// calling filter_pTmin(dummy=value)
//	// assumes that the declaration of a variable "float dummy;" is present
//	// at an appropriate point in the code of the filter.
//	//*************************************************
//
//	vector<vector<int> > outputclassmembers;
//	outputclassmembers.resize(0);
//	outputclassmembers = previousclassmemeber;
//
//	vector<string> assignments = getCSstrings(parameters);
//	map<string, float> mymap;
//	mymap = assignment2map(assignments);
//	// from here on the variables are initialized and valorized from the string provided as input
//	float etamin;
//	valorize(etamin, "etamin", mymap);
//
//	//---->loop over the particles of the input list !!!from the end!!!<------
//	for (Int_t jcount = previousclassmemeber.size() - 1; jcount > -1; --jcount) {
//
//		//do the checks that the user wants to do to select particles
//
//		vector<vector<int> > matrixtoadd;
//		matrixtoadd.resize(0);
//		matrixtoadd.push_back(previousclassmemeber[jcount]);
//
//		double eeta;
//		eeta = 0;
//		// call the observable that is needed
//		eeta = eta(matrixtoadd, vec_partvector, partvector);
//		// do the check on the computed value
//		if (eeta < etamin) {
//			outputclassmembers = stripelement(jcount, outputclassmembers);
//		}
//
//	}
//
//	return outputclassmembers;
//
//}

/*
 //*************************************************
 extern "C" vector<int> filter_autojet(vector<int> classmembersIN, vector<TRootLHEFParticle*> partvector, vector<TLorentzVector> vec_partvector) {
 //*************************************************
 vector<int> classmembersOUT;
 vector<int> jet;
 vector<int> notused;
 notused.resize(0);
 jet.resize(0);
 classmembersOUT.resize(0);

 notused = classmembersIN;

 long double Emax;
 Emax = 0;
 int mostenergeticparticle;
 int currentparticle;

 TLorentzVector jetfv;

 //---->loop over *particles*<------
 for (Int_t partincount = 0; partincount < classmembersIN.size(); ++partincount) {

 currentparticle = classmembersIN[partincount];

 if (vec_partvector[currentparticle].E() > Emax) {

 mostenergeticparticle = currentparticle;
 Emax = vec_partvector[currentparticle].E();

 }

 }

 //cout << "the most energetic particle is the PID " << partvector[mostenergeticparticle]->PID << " in row " << mostenergeticparticle << "with energy " << Emax << " GeV" << endl;

 moveparticle(mostenergeticparticle, notused, jet);

 //cout << "now the jet is made of ";
 //printvector(jet);
 //cout << "the particles left are ";
 //printvector(notused);


 for (int kappa = 0; kappa < 20 && kappa < notused.size(); ++kappa) {
 //---->loop over *particles*<------

 jetfv.SetPxPyPzE(0, 0, 0, 0);
 long double deltaRmin;
 deltaRmin = 100000000;
 int closestparticle;
 closestparticle = 10000000;

 for (int vect = 0; vect < jet.size(); ++vect) {

 jetfv = jetfv + vec_partvector[jet[vect]];

 }

 for (Int_t partincount = 0; partincount < notused.size(); ++partincount) {
 long double deltaR;
 currentparticle = notused[partincount];

 deltaR = vec_partvector[currentparticle].Vect().DeltaR(jetfv.Vect());
 if (deltaR < deltaRmin) {

 closestparticle = currentparticle;
 deltaRmin = deltaR;

 }

 }

 //cout<<" step: "<< kappa <<" deltaR: "<<
 cout << "{" << deltaRmin << "," << jetfv.E() << "},";

 moveparticle(closestparticle, notused, jet);

 }
 cout << endl;
 return jet;
 }
 */
////*************************************************
//extern "C" vector<vector<int> > with_pTmin(vector<vector<int> > inkplets,
//		vector<vector<int> > inclassmemeber,
//		vector<TLorentzVector> vec_partvector,
//		vector<TRootLHEFParticle*> partvector, string parameters, FUN OBS) {
//
//	//*************************************************
//	// the input string "parameters" must be a comma separated list without spaces
//	// of the form variable1=val1,variabe2=val2
//	// at run-time the use can specify an arbtrary number of variables
//	// but each of them must be explictly declared in the code of the filter he/she is using
//	// calling filter_pTmin(dummy=value)
//	// assumes that the declaration of a variable "float dummy;" is present
//	// at an appropriate point in the code of the filter.
//	//*************************************************
//
//
//	vector<string> assignments = getCSstrings(parameters);
//	map<string, float> mymap;
//	mymap = assignment2map(assignments);
//	// from here on the variables are initialized and valorized from the string provided as input
//	float pTmin;
//	valorize(pTmin, "pTmin", mymap);
//
//	vector<vector<int> > outputclassmembers;
//	outputclassmembers.resize(0);
//	outputclassmembers = inkplets;
//	// if this is run as a filter for kplets it is exactly what you want
//	// when run as filter on single particles it make sense anyways because in that case
//	// the two input lists are equal
//
//	//---->loop over the ordered list and take the requested elements
//	for (Int_t jcount = 0; jcount < inclassmemeber.size(); jcount++) {
//
//		//do the checks that the user wants to do to select particles
//		vector<vector<int> > matrixtoadd;
//		matrixtoadd.resize(0);
//		//		if (inkplets == inclassmemeber) {
//		//			//old case in which the filter is called in a "particle"-fashion
//		//			matrixtoadd.push_back(inclassmemeber[jcount]);
//		//		} else {
//		//			//			// if the kplets are subject to the filtering it is needed to translate the kplet in the event numbering scheme
//		//			//			 computeOBSonthekpletsimple(partvector, vec_partvector,
//		//			//					vector<vector<int> > allkplets,  inclassmembers, FUN OBS)
//		//			//			matrixtoadd.push_back(previousclassmemeber[jcount]);
//		//		}
//
//
//		//		matrixtoadd.push_back();
//
//		double pT;
//		pT = 0;
//		// call the observable that is needed
//		pT = pTofn(matrixtoadd, vec_partvector, partvector);
//
//		// do the check on the computed value
//		if (pT < pTmin) {
//			printmatrix(outputclassmembers);
//			outputclassmembers = stripelement(jcount, outputclassmembers);
//			cout << jcount << " particle has been removed because it had pT="
//					<< pT << endl;
//			printmatrix(outputclassmembers);
//		}
//
//	}
//
//	return outputclassmembers;
//
//}
////*************************************************
//extern "C" vector<vector<int> > take_Nth_in_pT(vector<vector<int> > kplets,
//		vector<vector<int> > previousclassmemeber,
//		vector<TLorentzVector> vec_partvector,
//		vector<TRootLHEFParticle*> partvector, string parameters, FUN OBS) {
//
//	//*************************************************
//	// the input string "parameters" must be a comma separated list without spaces
//	// of the form variable1=val1,variabe2=val2
//	// at run-time the use can specify an arbtrary number of variables
//	// but each of them must be explictly declared in the code of the filter he/she is using
//	// calling filter_pTmin(dummy=value)
//	// assumes that the declaration of a variable "float dummy;" is present
//	// at an appropriate point in the code of the filter.
//	//*************************************************
//
//	//*************************************************
//	// this filter eliminate all the particles but the one
//	// that is the n-th in the sorting according to pT
//	//*************************************************
//
//
//	vector<vector<int> > outputclassmembers;
//	outputclassmembers.resize(0);
//	//outputclassmembers = previousclassmemeber;
//
//	vector<string> assignments = getCSstrings(parameters);
//	map<string, float> mymap;
//	mymap = assignment2map(assignments);
//	// from here on the variables are initialized and valorized from the string provided as input
//	float nth;
//	valorize(nth, "nth", mymap);
//
//	//---->loop over the particles of the input list !!!from the end!!!<------
//	// build the necessary pairs of pair<pT,EID> and sort according to the long double pT
//
//	vector<long double> object_to_sort;
//
//	//cout << endl << "prima del sort" << endl;	printmatrix(previousclassmemeber);	cout << endl;
//
//	for (Int_t jcount = 0; jcount < previousclassmemeber.size(); jcount++) {
//
//		//do the checks that the user wants to do to select particles
//
//		vector<vector<int> > matrixtoadd;
//		matrixtoadd.resize(0);
//		matrixtoadd.push_back(previousclassmemeber[jcount]);
//
//		long double obs;
//		// call the observable that is needed
//		obs = pTofn(matrixtoadd, vec_partvector, partvector);
//		object_to_sort.push_back(obs);
//	}
//
//	//printvector(object_to_sort);	cout << endl;
//
//	vector<vector<int> > sorted_classmembers;
//	sorted_classmembers = sort2ndaccordingto1st(object_to_sort,
//			previousclassmemeber);
//
//	//cout << "dopo il sort" << endl;	printmatrix(sorted_classmembers);	cout << endl;
//
//	//	// print the sorted observables for a check
//	//	for (Int_t jcount = 0; jcount < sorted_classmembers.size(); jcount++) {
//	//
//	//		vector<vector<int> > matrixtoadd;
//	//		matrixtoadd.resize(0);
//	//		matrixtoadd.push_back(sorted_classmembers[jcount]);
//	//
//	//		long double obs;
//	//		// call the observable that is needed
//	//		cout << pTofn(matrixtoadd, vec_partvector, partvector) << " ";
//	//	}
//
//	// strip all the particle except the one that satisfies the ordering requirement
//	for (Int_t jcount = 0; jcount < sorted_classmembers.size(); jcount++) {
//		if (jcount + 1 == nth) {
//			outputclassmembers.push_back(sorted_classmembers[jcount]);
//		}
//
//	}
//
//	return outputclassmembers;
//
//}
#endif
